function [beta,se,p,r2,F,control_mean,N,C,T] = ipwreg(y,treat,X,weight,stations,tail);
    %%Estimation Sample%%
    [N,k] = size(X);
    T = sum(treat);
    C = sum(1-treat);
    
    weight_sqrt = sqrt(weight);
    
    control_mean = sum((1-treat).*weight.*y)/sum((1-treat).*weight);
    
    Xw = [weight_sqrt.*treat repmat(weight_sqrt,1,k).*X];
    yw = weight_sqrt.*y;
    [beta,bint,e,rint,stats] = regress(yw,Xw);
    
    %beta = inv(Xw'*Xw)*Xw'*yw;
    %e = yw-Xw*beta;
    
    dep_matrix = depmat(N,stations);
        
    omega = (e*e').*dep_matrix;
    V = inv(Xw'*Xw)*(Xw'*omega*Xw)*inv(Xw'*Xw);
    se = sqrt(diag(V));
    
    sst = sum((yw-mean(yw)).^2);
    ssm = sum((Xw*beta-mean(yw)).^2);
    sse = sum(e.^2);
    
    %r2 = ssm/sst;
    r2 = stats(1);
    F = (ssm/k)/(sse/(57-k-1));
    
    p = zeros(k+1,1);
    
    for i = 1:k+1
        if tail == 1
            p(i) = 1 - cdf('norm',abs(beta(i)/se(i)),0,1);
        elseif tail == 2 
            p(i) = 2*(1 - cdf('norm',abs(beta(i)/se(i)),0,1));
        end;
    end;

    
    
        
    
    
    
    
    
    